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Abstract. A new dynamical clustering algorithm for the identification of modules in complex 
networks has been recently introduced fll. In this paper we present a modified version of this 
algorithm based on a system of chaotic Rossler oscillators and we test its sensitivity on real and 
. computer generated networks with a well known modular structure. 
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INTRODUCTION 



An important property common to many networks is the presence of modules or com- 
munity structures, that can be roughly defined as subsets of network nodes within which 
the network connections are dense, but between which they are sparser. Since nodes 
belonging to tightly connected clusters of nodes are most likely to have other properties 
in common, the detection of these structures in complex networks is potentially very 
useful. 

We have recently proposed a dynamical clustering (DC) method for the modules iden- 
tification based on the properties of a dynamical system associated to the graph [Q]]. 
Such a method combines topological and dynamical information in order to devise an 
algorithm that is able to identify the modular structure of a graph with a precision and 
a computational cost (0(N 2 )) competitive with the best techniques based solely on the 
topology. The method is based upon the well-known phenomenon of synchronization 
clusters of non identical phase oscillators [0], each one associated to a node, and inter- 
acting through the edges of the graph. Clusters of synchronized oscillators represent an 
intermediate regime between global phase locking and full absence of synchronization, 
thus implying a division of the whole graphs into groups of elements which oscillate at 
the same (average) frequency. The key idea is that, starting from a fully synchronized 
state of the network, a dynamical change in the weights of the interactions, retaining 
information of the original betweenness distribution, yields a progressive hierarchical 
clustering that fully detects modular structures. 

In this paper we implement our algorithm on a systems of N identical (three- 
dimensional) chaotic Rossler oscillators and we test the precision (sensitivity) obtained 



by the algorithm on real and computer generated networks with well known modular 
structures. 

DYNAMICS OF A WEIGHTED NETWORK OF ROSSLER 

OSCILLATORS 

The dynamics of a network of N coupled identical oscillators is described by: 

N 

±i = F(x j )-CT£GyH[x I — Xy], i=l,...,N, (1) 
;'=i 

where F(x) governs the dynamics of each individual oscillator, H(x) is a linear vectorial 
function, o is the overall coupling strength and the G = Gjj is the coupling matrix. The 
rows of matrix G have zero sum and this ensures that the completely synchronized state 
{xi(t) = s(?),Vi | s = F(s)} is a solution of Eq. ([]])• By means of the so called Master 
Stability Function approach, it is possible to study the conditions under which such a 
state is stable [3], i.e. the propensity for synchronization (PFS) of a given network. 
In Ref.[3] it was shown that an enhancement in the PFS can be achieved by exploiting 
the information contained in the overall topology of the network. This can be done 
through an opportune choice of the coupling matrix G that makes use of the load concept 
and by scaling the coupling strenght a in Eq. © to the load of each link. The load of 
the link connecting nodes i and j is quantified by the so called edge betweenness, i.e. the 
fraction of shortest paths that are making use of that link. By means of this weighting 
procedure, that clearly reflects the network structure at a global scale, Eq. (OQ) reads: 

** = FW-rnaE?^- 1 /] j=i,...,iv, (2) 

where a is a real tunable parameter, and Kj is the set of neighbors of node i th . For a 
given dynamical system F(x ; ), for a given value of o and for a given network topology 
it is possible to find, by means of the Master Stability Function approach what is the 
value of ttbest providing the best PFS of the system. In practice, it is more convenient to 
put a = and to find a value of the coupling parameter o which would ensure a fully 
synchronized state for the oscillators network. 

In Ref.[l] we used as dynamical system the so-called Opinion Changing Rate (OCR) 
model with an Heigselmann-Krause dynamics qO, and we showed that, starting 
from a perfectly synchronized state for a = 0, if a is let to decrease from to — °°, the 
links with the higher load will be weighted less and less with respect to the other links, 
thus inducing a progressive desynchronization of the system in clusters of frequencies 
(dynamical clustering) corresponding to different modules, or communities, of a given 
network. Here we apply our analysis to real or trial networks using a system of chaotic 
Rossler oscillators and we study again the dynamical clustering process as a function of 
decreasing values of a, identifying a likely community subdivision of the networks by 
looking to local or global maxima of the modularity Q[l]. The latter simply quantifies 
the degree of correlation between the probability of having an edge joining two sites and 



FIGURE 1. The Karate Club network , with the two Zachary's communities identified by circles and 
squares. 



the fact that the sites belong to the same community [6], thus in general it makes sense 
to look for large values of Q. In fact we get Q = if we consider the whole network as 
a single community or if we consider a completely random network. On the other hand, 
for networks with an appreciable subdivision in classes, Q usually falls in the range 
between 0.2 and 0.7. 

The dynamics of a system of /V identical (three-dimensional) chaotic Rossler oscillators, 
defined over the nodes of a given network, is ruled by Eq.©, with x ; - = F(x,) = 

[— (Dyi — Zi,(OXi + 0.1 65y/, 0.2+ Zifa — 10)] and H(x) = [x,0,0] (thus the coupling acts 
only on the x variable). In other words we have the following equations of motion: 

it = -coyi-n - ° £ If} (xi -Xj) 

L l K- ''/ / K, 

ji = (OXi + 0A65yi (3) 

Zi = 0.2 + 10) *=1,...,/V, 

Here ft) is a common natural frequency associated at each oscillator that, without loss 
of generality, we put equal to 1 .0. The load matrix kj (the matrix of the edge between- 
nesses) is calculated once forever for the chosen network with a computational cost of 
O(KN), K being the total number of links. 

In order to evaluate the degree of synchronization of the Rossler system () one has to cal- 
culate the order parameter *F = (^|£jLi e^ l ^\) t , where 4> ; (?) = arctan[j^} indicates 
the istantaneous phase of the ?'-th oscillator and (...)/ stays for a time average. If all the 
oscillators rotate independently, no clusters exist and we have *F ~ -^U. On the contrary, 
if their motions are synchronized in phase, only one cluster exists and we obtain *F ~ 1 . 
Once a network is fixed, the first task is to find the value of the coupling parameter o 
providing a fully synchronized starting state for the Rossler oscillators at a = (i.e. at 
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FIGURE 2. Asymptotic Rtissler order parameter (at a = 0) versus the coupling a for the Karate Club 
network. See text for further details. 



t = 0). Then, one can let a to decrease in time and study the dynamical clustering pro- 
cess acting on the istantaneous phases <!>;(?) 's of the oscillators. Notice that these phases 
play here the same role played by the istantaneous frequencies in the OCR model [Q]]: 
in this case we call "cluster" a group of contiguous phases in the 4>'s interval (usually 
[—3,3]) separated by a distance of more than 0.02 units. For each value of a a different 
configuration of clusters (corresponding to a given network structure) will appear and 
one has to calculate the corresponding modularity and select the configuration with the 
best modularity score. 

In the following we will show in detail this process by putting the Rossler system over 
different real and trial networks. 



Zachary's karate club 

As first example, in order to test our algorithm for finding community struc- 
tures, we consider a real network, the well-known Karate Club network analyzed 
by Zachary fl7D. It consists of N = 34 individuals (nodes),- whose mutual friend- 
ship relations (expressed by K = 78 edges) have been carefully investigated over 
a period of two years. Due to contrasts between a teacher and the administrator 
of the club, the club splitted into two smaller communities. The corresponding 
network is presented in Fig.l, where squares and circles label the members of 
the two groups. The 'circles' community has 18 elements (corresponding to nodes 
9,10,15,16,19,21,23,24,25,26,27,28,29,30,31,32,33,34) while the 'squares' com- 
munity has 16 elements, (nodes 5,6,7,11,17,1,2,3,4,8, 12,13,14,18,20,22). The 
a- priori modularity of such a configuration results to be equal to <2z ~ 0.37. 
Firstly, in Fig. 2 we plot the behavior of the asymptotic Rossler order parameter (aver- 
aged over 10 events) versus the coupling strenght o for the Karate Club network with 
a = 0. It results that, above o ~ 0.6, the system lies in the fully synchronized phase, 
thus in the following we reasonably set a = 1.3. Using such a value in the equations of 
motion (), we can now integrate them numerically and study the dynamical clustering 
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FIGURE 3. Typical run for the Karate Club Network: time-evolution of both the Rossler's phases 
configurations (top panel) and the corresponding modularity (bottom panel) as a function of a(t). See 
text for further details. 



TABLE 1. Clusters configuration with the best modularity score 
Qbest = 0.42 for the Karate Club Network (at a best ~ -1.76). See 
text for further details. 

cluster 1 (11 nodes) 1,2,3,4,8,12,13,14,18,20,22 

cluster 2 (12 nodes) 9,10,15,16,19,21,23,27,29,30,31,32,33,34 

cluster 3 (4 nodes) 24,25,26,28 

cluster 4 (2 nodes) 27,30 

cluster 5 (5 nodes) 5,6,7,11,17 



desynchronization process acting on the oscillators' phases. The task of our algorithm 
is to extract the best modular structure of the network by using only the information 
expressed by the edge betweenness of its links, which are calculated only once for the 
given network. 

The system starts in a perfectly synchronized state for a start = 0, i.e. jc,- (0) = y;(0) = 
Zi(0) = Vz, thus 4> ; (0) = Vz and only one cluster exists (choosing a different starting 
value for Xi,yi,Zi does not change significantly the conclusions shown in the following); 
then we decrease a during a single run as a function of time with a rate of 2 time steps 
and a given constant decrement da = 0.0008. Simultaneously, for each value of a(t) 
and for each configuration of clusters, the modularity Q is calculated. 
In the top panel of Fig. 3 we show the phases time evolution as a function of a(t) in 
a typical run with a going from to —3. In this figure 34 points are plotted for each 
value of a(t), corresponding to the istantaneous phases of the oscillators (from which 
the average istantaneous phase of the system has been subtracted in order to have a 
symmetric plot). Correspondingly, in the bottom panel of Fig. 3, the modularity Q is 
also plotted as a function of a(t). Even if the system strongly oscillates during the 



desynchronization process (at variance with the much more stable behavior of the OCR 
system in [1]), clusters configurations (i.e. community structures of the underlying 
network) with very large values of modularity appear. 

In Table 1 we report the detailed clusters configuration corresponding to the maximum 
value of modularity, Qb es t = 0.42, found for a^est ~ —1 .76 (see the arrow in the bottom 
panel). One can see that five clusters have been found: the first and the last ones, 
respectively made of 11 and 5 nodes (oscillators), if considered together, correspond 
to the first community of 16 nodes (squares in the left panel of Fig.l) observed by 
Zachary, while the sum of the remaining three clusters corresponds to tightly connected 
subgroups of the second community of 18 nodes (circles in the left panel of Fig.l). On 
the other hand, for several values of a in Fig. 3 (e.g. for a = —0.538 or a = —1.075) 
the clusters configuration corresponding to the two Zachary communities has been also 
recovered, but - as previously seen -, being its modularity smaller than 0.42 (in fact 
Qz ~ 0.37), the algorithm favours the five clusters configuration shown in Table 1 (in 
other words, the Zachary community subdivision corresponds only to a local maximum 
of modularity for this particular network and not to a global one). 
In conclusion, at least in so far as it concerns the Zachary Club network, the dynamical 
clustering algorithm based on the Rossler system seems to work very well, even if 
compared with the analogous results presented in JH]. In fact it results to be at the same 
time very fast (it extracts the best clusters configuration of the network versus a in a 
single run) and very sensitive (it is able to recover community structures of the Zachary 
network with a great modularity). In the next section we will apply the same algorithm 
to another real network, the Chesapeake Bay food web. 



Chesapeake Bay Food Web 

Another classical benchmark for the community identification algorithms is a food 
web of marine organisms living in the Chesapeake Bay, situated on the Atlantic coast 
of the United States. This ecosystem was originally studied by Baird and Ulanowicz 
dsl], who carefully investigated the predatory interactions between the most important 
taxa (species or groups of species) and constructed a network of 33 vertex and 71 links. 
We will consider here (as usually done in many papers lESGll]) its non-directed and 
non-valued version and calculate first (only once before each simulation using the same 
network) the edge betweenness of each link. 

In Fig.4 we show the result of a typical event obtained in a simulation performed 
with the same procedure described in the previous section and for a value of the in- 
teraction strenght o = 1 (such that the Rossler's system would lie in its synchro- 
nized phase for a — 0). The system starts in a fully synchronized state (we set again 
jcj(0) = y;(0) = n(0) = V/) at a start = and evolves through decreasing values of a(t) 
(with a decrement 8a = 0.0008), up to the value a en d = —2. The clusters evolution 
and the corresponding modularity Q(t) are plotted as a function of a and, again, look 
very oscillating in time. The detailed configuration with the highest modularity peak 
(see the arrow in the bottom panel) is reported in Table 2. It consists of 6 clusters with a 
Qmax = 0.43, obtained for OfAvr = —1-62, and it is quite consistent with the correspond- 
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FIGURE 4. Typical run for the Chesapeake Bay food web: time-evolution of both the Rossler's phases 
configurations (top panel) and the corresponding modularity (bottom panel) as a function of a(t). 



TABLE 2. Clusters configuration with the best modularity 
score Qbest = 0.43 for the Chesapeake Bay food web (at at, est ~ 
-1.62) 



cluster 1 (10 nodes) 


3,14,15,16,18,25,26,27,28,29 


cluster 2 (3 nodes) 


4,17,19 


cluster 3(1 nodes) 


30 


cluster 4 (3 nodes) 


22,31,32 


cluster 5(14 nodes) 


1,2,7,8,9,10,11,12,13,20,21,23,24,33 


cluster 6 (2 nodes) 


5,6 



ing results of [Ql |9|, Oil]), where a main separation in two large communities had been 
found, according to the distinction between pelagic organisms, which live near the sur- 
face or at middle depths (clusters n.1,2 and 4), and benthic organisms, which live near 
the bottom (clusters n.3,5 and 6). 

Such a result corroborates the good performance of the Rossler algorithm in the identifi- 
cation of community structures in real networks. The next step will be to test this method 
on "ad hoc" trial networks with a well known fixed community structure OLSl, in order 
to explore in deeper detail its sensitivity. 

Sensitivity test for ad hoc trial networks 

Typical trial networks are generated with N = 12& nodes and split into four commu- 
nities containing 32 nodes each. Pairs of nodes belonging to the same community are 
linked with probability /?,„, whereas pairs belonging to different communities are joined 
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FIGURE 5. Sensitivity curves of the dynamical clustering algorithm with Rossler system: the fraction 
of correctly classified nodes in a set of N=128 trial networks is reported as a function of z out for several 
increasing values of the decrement 8a. 



with probability p out . The value of p out is taken so that the average number of links a 
node has to members of any other community, i.e. z ou t, can be controlled. While p out 
(and therefore Zout) is varied freely, the value of p,„ is chosen to keep the total average 
node degree k constant, and set to 16. As z ut is increased from 2 (very well defined 
structures) to 8 (bad defined structures), the communities become more and more dif- 
fuse and harder to detect. Since the "real" community structure is well known in this 
case, it is possible to measure the number of nodes correctly classified by our method of 
community identification (see for example yj]). 

We apply our algorithm to 10 different sets of trial networks. For each network, as done 
in the previous sections, after having calculated the load matrix Uj, we integrate numeri- 
cally the equations of motion (). Every time we start from a perfectly synchronized state 
for OCstart = and we analyze the desynchronization process when a(t) decreases in time 
(with a given constant decrement 8a) during a single run. In Fig. 5 we plot the number 
of correctly classified nodes as a function of z ut, averaged over the set of 10 networks 
and for three increasing values of 8a. The error bars (standard deviations) for each point 
is also reported. As one can see, the smaller is 8a, the better is the result. However, in 
any case the sensitivity abruptly falls above z ut = 5, staying around the 40% of correctly 
identified nodes up to Zout = 8. Such a performance of the Rossler system is surely worse 
than that of the OCR — HK system shown in yj], but in any case it confirms the possibil- 
ity of extending the dynamical clustering algorithm presented in yj] to other dynamical 
systems (also three-dimensional, like in this case) with quite good results. 

CONCLUSION 

Summarizing, even if the global sensitivity of the dynamical clustering (DC) algorithm 
for the a system of chaotic Rossler oscillators seems to be not competitive with respect 



to those of other methods (see Refs.O], O, [12]), on the other hand the results here 
presented confirm that (i) the DC algorithm is robust with respect to the change of the 
adopted dynamical system; (ii) the system of chaotic Rossler oscillators works quite well 
if applied to the real networks considered; (iii) the DC algorithm is quite fast and needs 
only few runs of integration for each network (after the calculation of the corresponding 
load matrix Uj) imAWm- 
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